
/*
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */


#include <stdio.h>
#include <limits.h>

#include "libmooigraph.h"
#include "lmg.h"

#define MAX_LEVELS 1024

typedef struct
{
  int redirection;
  float barycenter;
  struct lmgnode *n;
} redirection_element;

typedef struct
{
  int rows;
  int columns;
  int nedges;
  redirection_element *logical_to_physical_row_number;
  redirection_element *logical_to_physical_col_number;
  /* bitset set to 1 if there is an edge between row and column */
  unsigned char *elembits;
} matrix;


static void privateUpdateNodePositions (int debug);
static void privateBuildAndInitializeMatrices (void);
static void matrixArrayDownUpIterate (matrix * m[], int size,
				      int max_iterations);
static unsigned int matrixArrayCrossings (matrix * array[], int size);
static void matrixArrayDownUp (matrix * m[], int size);
static unsigned int current_edge_crossings;
static matrix *matrices[MAX_LEVELS];
static matrix *matrixCreate (int r, int c);
static void matrixAtRowPutNode (matrix * m, int r, struct lmgnode *n1);
static void matrixAtColumnPutNode (matrix * m, int c, struct lmgnode *n1);
static void matrixAtPut (matrix * m, int r, int c, int element);
static struct lmgnode *matrixNodeAtRow (matrix * m, int r);
static struct lmgnode *matrixNodeAtColumn (matrix * m, int c);
static int nswap;


void
lmg_barycenter (void)
{
  int crossingIterations;
  nswap = 0;
  crossingIterations = 10;

  if (lmgdebug)
    {
      printf ("%s(): start\n", __FUNCTION__);
    }

  /* build and initialize a matrix representation of the selected
     nodes and their connecting edges */
  privateBuildAndInitializeMatrices ();

  if (lmgdebug)
    {
      printf ("%s(): iterate\n", __FUNCTION__);
    }

  matrixArrayDownUpIterate (matrices, lmgnlevels, crossingIterations);

  /* update the nodes' position variable to their position in the
     matrices after the permutation */
  privateUpdateNodePositions (lmgdebug);

  if (lmgdebug)
    {
      printf ("%s(): finish\n", __FUNCTION__);
    }

  return;
}

void
lmg_barycenterclear (void)
{
  matrix *res;
  int i;
  if (matrices == NULL)
    {
      return;
    }

  for (i = 0; i < lmgnlevels; i++)
    {
      res = matrices[i];
      if (res)
	{
	  lmg_free (res->logical_to_physical_row_number);
	  lmg_free (res->logical_to_physical_col_number);
	  lmg_free (res->elembits);
	  lmg_free (res);
	}
      matrices[i] = NULL;
    }

  return;
}

unsigned int
lmg_crossings (void)
{
  return (current_edge_crossings);
}

static void
privateBuildAndInitializeMatrices (void)
{
  int i;
  struct lmgnodel *ptr;
  struct lmgnode *n1;
  struct lmgnode *n2;
  struct lmgedgel *el;
  for (i = 0; i < lmgnlevels; i++)
    {
      /* create level matrix */
      matrices[i] = matrixCreate (lmgnal[i], lmgnal[i + 1]);
    }

  for (i = 0; i < lmgnlevels; i++)
    {
      ptr = lmglevels[i];
      while (ptr)
	{
	  n1 = ptr->node;
	  matrixAtRowPutNode (matrices[i], n1->tree_layout.position, n1);
	  ptr = ptr->next;
	}
    }

  for (i = 0; i < lmgnlevels; i++)
    {
      ptr = lmglevels[i + 1];
      while (ptr)
	{
	  n1 = ptr->node;
	  matrixAtColumnPutNode (matrices[i], n1->tree_layout.position, n1);
	  ptr = ptr->next;
	}
    }

  el = lmgedgelist;
  while (el)
    {
      n1 = el->edge->from;
      n2 = el->edge->to;
      i = n1->tree_layout.level;
      if (n1->tree_layout.level + 1 == n2->tree_layout.level)
	{
	  matrixAtPut (matrices[i], n1->tree_layout.position,
		       n2->tree_layout.position, 1);
	}
      else
	{
	  printf
	    ("%s(): skipping edge between '%s' level=%d and '%s' level=%d\n",
	     __FUNCTION__, n1->name, n1->tree_layout.level, n2->name,
	     n2->tree_layout.level);
	}
      el = el->next;
    }

  return;
}


static void
privateUpdateNodePositions (int debug)
{
  int i;
  int j;
  struct lmgnode *n1;
  int number_of_nodes_at_level;
  const char *s = NULL;
  for (i = 0; i < lmgnlevels; i++)
    {
      number_of_nodes_at_level = lmg_nal (i);
      for (j = 0; j < number_of_nodes_at_level; j++)
	{
	  n1 = matrixNodeAtRow (matrices[i], j);
	  if (n1)
	    {
	      if (debug)
		{
		  if (n1->tree_layout.position != j)
		    {
		      if (n1->name)
			{
			  s = n1->name;
			}
		      else
			{
			  s = "dummy";
			}
		      printf ("%s(): moving node %s from position %d to %d\n",
			      __FUNCTION__, s, n1->tree_layout.position, j);
		    }
		}
	      n1->tree_layout.position = j;
	    }
	}
    }
  i = lmgnlevels;
  number_of_nodes_at_level = lmg_nal (i);
  for (j = 0; j < number_of_nodes_at_level; j++)
    {
      n1 = matrixNodeAtColumn (matrices[i - 1], j);
      if (n1)
	{
	  if (debug)
	    {
	      if (n1->tree_layout.position != j)
		{
		  if (n1->name)
		    {
		      s = n1->name;
		    }
		  else
		    {
		      s = "dummy";
		    }
		  printf ("%s(): moving node %s from position %d to %d\n",
			  __FUNCTION__, s, n1->tree_layout.position, j);
		}
	    }

	  n1->tree_layout.position = j;
	}
    }
  return;
}

/* iterate "down" and "up" procedures to reduce the number of edge crossings */
static void
matrixArrayDownUpIterate (matrix * m[], int size, int max_iterations)
{
  int i = 0;
  unsigned int oldcross = 0;
  unsigned int startedcross = 0;
  if (lmgdebug)
    {
      printf ("%s(): start size=%d\n", __FUNCTION__, size);
    }
  /* set to -1 to indicate unknown */
  current_edge_crossings = UINT_MAX;
  oldcross = 0;
  if (lmgdebug)
    {
      printf ("%s(): starting barycenter with %d edge line crossings\n",
	      __FUNCTION__, current_edge_crossings);
    }
  for (i = 0; i < max_iterations; i++)
    {
      if (current_edge_crossings == 0)
	{
	  if (lmgdebug)
	    {
	      printf ("%s(): After iteration %d, crossings: %d\n",
		      __FUNCTION__, i, current_edge_crossings);
	    }
	  if (lmgdebug)
	    {
	      printf ("%s(): finish barycenter with %d edge line crossings\n",
		      __FUNCTION__, current_edge_crossings);
	    }
	  return;
	}
      matrixArrayDownUp (m, size);
      current_edge_crossings = matrixArrayCrossings (m, size);
      if (startedcross == 0)
	{
	  startedcross = current_edge_crossings;
	}
      if (lmgdebug)
	{
	  printf ("%s(): After iteration %d, crossings: %d nswap=%d\n",
		  __FUNCTION__, i, current_edge_crossings, nswap);
	}
      if (oldcross == current_edge_crossings)
	{
	  break;
	}
      oldcross = current_edge_crossings;
    }
  if (lmgdebug)
    {
      printf
	("%s(): finish barycenter with %d edge line crossings and started with %d edge line crossings\n",
	 __FUNCTION__, current_edge_crossings, startedcross);
    }
  return;
}

/* edge bits */
#define BIT_SET(  p, n) (p[(n)/CHAR_BIT] |=  (0x80>>((n)%CHAR_BIT)))
#define BIT_CLEAR(p, n) (p[(n)/CHAR_BIT ] &= ~(0x80>>((n)%CHAR_BIT )))
#define BIT_ISSET(p, n) (p[(n)/CHAR_BIT] &   (0x80>>((n)%CHAR_BIT)))

/* return the element at (r,c) */
static inline int
matrixAt (matrix * m, int r, int c)
{
  int i = 0;
  int ret = 0;
  i = m->logical_to_physical_row_number[r].redirection * m->columns
    + m->logical_to_physical_col_number[c].redirection;
  if (BIT_ISSET (m->elembits, i))
    {
      ret = 1;
    }
  else
    {
      ret = 0;
    }
  return (ret);
}

/* assign 'element' to (r,c) set to 1 */
static void
matrixAtPut (matrix * m, int r, int c, int element)
{
  int i = 0;

  i = m->logical_to_physical_row_number[r].redirection * m->columns
    + m->logical_to_physical_col_number[c].redirection;

  if (element)
    {
      BIT_SET (m->elembits, i);
    }
  else
    {
      BIT_CLEAR (m->elembits, i);
    }

  /* added an edge to the matrix */
  m->nedges = m->nedges + 1;

  return;
}

/* associate the given node with the matrix row r */
static void
matrixAtRowPutNode (matrix * m, int r, struct lmgnode *n1)
{
  m->logical_to_physical_row_number[r].n = n1;
  return;
}

/* associate the given node with the matrix column c */
static void
matrixAtColumnPutNode (matrix * m, int c, struct lmgnode *n1)
{
  m->logical_to_physical_col_number[c].n = n1;
  return;
}

/* return the node associated with row r */
static struct lmgnode *
matrixNodeAtRow (matrix * m, int r)
{
  return m->logical_to_physical_row_number[r].n;
}

/* return the node associated with column c */
static struct lmgnode *
matrixNodeAtColumn (matrix * m, int c)
{
  return m->logical_to_physical_col_number[c].n;
}

/* deallocate matrix */
static void
matrixDestroy (matrix * m)
{
  lmg_free (m->logical_to_physical_row_number);
  lmg_free (m->logical_to_physical_col_number);
  lmg_free (m->elembits);
  lmg_free (m);
  return;
}

/* create a new matrix identical to m. */
static inline matrix *
matrixCopy (matrix * m)
{
  matrix *res;
  int i;
  int r;
  int c;
  int nbyte = 0;

  r = m->rows;
  c = m->columns;

  res = matrixCreate (r, c);

  for (i = 0; i < r; i++)
    {
      res->logical_to_physical_row_number[i] =
	m->logical_to_physical_row_number[i];
    }
  for (i = 0; i < c; i++)
    {
      res->logical_to_physical_col_number[i] =
	m->logical_to_physical_col_number[i];
    }
  nbyte = ((r * c) + CHAR_BIT) / CHAR_BIT;
  for (i = 0; i < nbyte; i++)
    {
      res->elembits[i] = m->elembits[i];
    }
  return res;
}

/* calculate the row barycenter of r */
static float
matrixRowBarycenter (matrix * m, int r)
{
  int l;
  int q;
  int sum1;
  int sum2;

  q = m->columns;
  sum1 = sum2 = 0;
  for (l = 0; l < q; l++)
    {
      sum1 += (l + 1) * matrixAt (m, r, l);
      sum2 += matrixAt (m, r, l);
    }
  if (sum1 == 0 || sum2 == 0)
    {
      return 0.0;
    }
  return (float) sum1 / (float) sum2;
}

/* assign barycenter to row */
static void
matrixAtRowPutBarycenter (matrix * m, int r, float b)
{
  int i;
  i = m->logical_to_physical_row_number[r].redirection;
  m->logical_to_physical_row_number[i].barycenter = b;
  return;
}

/* assign barycenter to column */
static void
matrixAtColumnPutBarycenter (matrix * m, int c, float b)
{
  int i;
  i = m->logical_to_physical_col_number[c].redirection;
  m->logical_to_physical_col_number[i].barycenter = b;
  return;
}

/* calculate the column barycenter of c */
static float
matrixColumnBarycenter (matrix * m, int c)
{
  int k;
  int p;
  int sum1;
  int sum2;

  p = m->rows;
  sum1 = 0;
  sum2 = 0;
  for (k = 0; k < p; k++)
    {
      sum1 = sum1 + (k + 1) * matrixAt (m, k, c);
      sum2 = sum2 + matrixAt (m, k, c);
    }
  if (sum1 == 0 || sum2 == 0)
    {
      return 0.0;
    }
  return ((float) sum1 / (float) sum2);
}

/* calculate row barycenters of all rows in m */
static void
matrixComputeRowBarycenters (matrix * m)
{
  int r;
  for (r = 0; r < m->rows; r++)
    {
      matrixAtRowPutBarycenter (m, r, matrixRowBarycenter (m, r));
    }
  return;
}

/* calculate column barycenters of all columns in m */
static void
matrixComputeColumnBarycenters (matrix * m)
{
  int c;
  for (c = 0; c < m->columns; c++)
    {
      matrixAtColumnPutBarycenter (m, c, matrixColumnBarycenter (m, c));
    }
  return;
}

/* calculate row and column barycenters of m
 * static void
 * matrixComputeBarycenters (matrix * m)
 * {
 *   matrixComputeRowBarycenters (m);
 *   matrixComputeColumnBarycenters (m);
 *   return;
 * }
 */

/* swap elements i and j in the redirection vector 'array[]' */
static void
swapElements (redirection_element array[], int i, int j)
{
  redirection_element tmp;
  nswap = nswap + 1;
  tmp = array[i];
  array[i] = array[j];
  array[j] = tmp;
  return;
}

/* bubble sort the redirection vector in order of ascending barycenters */
static void
sortRedirectionVectorOnBarycenter (redirection_element array[],
				   int array_size)
{
  int i;
  int j;

  for (i = 1; i < array_size; i++)
    {
      for (j = i; j > 0 && array[j - 1].barycenter > array[j].barycenter; j--)
	{
	  swapElements (array, j - 1, j);
	}
    }
  return;
}

/* bubble sort the redirection vector in order of ascending barycenters.
   Keep the position of elements in the vector 'follow_array[]' relative to
   the vector being sorted ('array[]') constant */
static void
sortRedirectionVectorOnBarycenterFollow (redirection_element array[],
					 redirection_element follow_array[],
					 int array_size)
{
  int i;
  int j;

  for (i = 1; i < array_size; i++)
    {
      for (j = i; j > 0 && array[j - 1].barycenter > array[j].barycenter; j--)
	{
	  swapElements (array, j - 1, j);
	  swapElements (follow_array, j - 1, j);
	}
    }
  return;
}

/* sort the rows of m in order of ascending barycenters. */
static void
matrixSortRowsOnBarycenter (matrix * m)
{
  sortRedirectionVectorOnBarycenter (m->logical_to_physical_row_number,
				     m->rows);
  return;
}

/* sort the columns of m in order of ascending barycenters. */
static void
matrixSortColumnsOnBarycenter (matrix * m)
{
  sortRedirectionVectorOnBarycenter (m->logical_to_physical_col_number,
				     m->columns);
  return;
}

/* sort the rows of m in order of ascending barycenters. 
   Keep the order of the columns in the matrix 'follow' relative to
   the order of the rows in 'm' constant */
static void
matrixSortRowsOnBarycenterFollow (matrix * m, matrix * follow)
{
  if (m->rows != follow->columns)
    {
      printf ("Dimension mismatch in matrixSortRowsOnBarycenterFollow()\n");
      return;
    }
  sortRedirectionVectorOnBarycenterFollow (m->logical_to_physical_row_number,
					   follow->logical_to_physical_col_number,
					   m->rows);
  return;
}

/* sort the columns of m in order of ascending barycenters. 
   Keep the order of the rows in the matrix 'follow' relative to
   the order of the columns in 'm' constant */
static void
matrixSortColumnsOnBarycenterFollow (matrix * m, matrix * follow)
{
  if (m->columns != follow->rows)
    {
      printf
	("Dimension mismatch in matrixSortColumnsOnBarycenterFollow()\n");
      return;
    }
  sortRedirectionVectorOnBarycenterFollow (m->logical_to_physical_col_number,
					   follow->logical_to_physical_row_number,
					   m->columns);
  return;
}

/* "down" procedure for reducing the number of edge crossings */
static inline void
matrixArrayDown (matrix * m[], int size)
{
  int i;
  unsigned int c1 = 0;
  unsigned int c2 = 0;
  matrix *m1[MAX_LEVELS];

  if (size >= MAX_LEVELS)
    {
      return;
    }
  /* save a "backup" of the matrices, in case the "down" procedure 
     makes the result _worse_ (more edge crossings) rather than better */
  for (i = 0; i < size; i++)
    {
      m1[i] = matrixCopy (m[i]);
    }

  c1 = current_edge_crossings;
  if (lmgdebug && 0)
    {
      printf ("%s(): In down, before: %d\n", __FUNCTION__, c1);
    }
  for (i = 0; i < size - 1; i++)
    {
      matrixComputeColumnBarycenters (m[i]);
      matrixSortColumnsOnBarycenterFollow (m[i], m[i + 1]);
    }
  matrixComputeColumnBarycenters (m[i]);
  matrixSortColumnsOnBarycenter (m[i]);

  c2 = matrixArrayCrossings (m, size);
  if (lmgdebug && 0)
    {
      printf ("%s(): In down, crossings after: %d\n", __FUNCTION__, c2);
    }
  if (c2 <= c1)
    {				/* if the result has been improved */
      current_edge_crossings = c2;
      for (i = 0; i < size; i++)	/* destroy backup */
	{
	  matrixDestroy (m1[i]);
	}
    }
  else
    {
      for (i = 0; i < size; i++)
	{			/* otherwise destroy the original */
	  matrixDestroy (m[i]);	/*   and use the backup           */
	  m[i] = m1[i];
	}
    }
  return;
}

/* "up" procedure for reducing the number of edge crossings */
static inline void
matrixArrayUp (matrix * m[], int size)
{
  int i;
  unsigned int c1 = 0;
  unsigned int c2 = 0;
  matrix *m1[MAX_LEVELS];

  if (size >= MAX_LEVELS)
    {
      return;
    }
  /* save a "backup" of the matrices, in case the "up" procedure 
     makes the result _worse_ (more edge crossings) rather than better */
  for (i = 0; i < size; i++)
    {
      m1[i] = matrixCopy (m[i]);
    }

  c1 = current_edge_crossings;
  if (0)
    {
      printf ("In up, before: %d\n", c1);
    }
  for (i = size - 1; i > 0; i--)
    {
      matrixComputeRowBarycenters (m[i]);
      matrixSortRowsOnBarycenterFollow (m[i], m[i - 1]);
    }

  matrixComputeRowBarycenters (m[0]);
  matrixSortRowsOnBarycenter (m[0]);

  c2 = matrixArrayCrossings (m, size);
  if (lmgdebug && 0)
    {
      printf ("%s(): In up, after: %d\n", __FUNCTION__, c2);
    }
  if (c2 <= c1)
    {				/* if the result has been improved */
      current_edge_crossings = c2;
      for (i = 0; i < size; i++)	/* destroy the backup */
	{
	  matrixDestroy (m1[i]);
	}
    }
  else
    {
      for (i = 0; i < size; i++)
	{			/* otherwise destroy the original */
	  matrixDestroy (m[i]);	/*   and use the backup           */
	  m[i] = m1[i];
	}
    }
  return;
}


static void
matrixArrayDownUp (matrix * m[], int size)
{
  matrixArrayDown (m, size);
  if (current_edge_crossings == 0)
    {
      return;
    }
  matrixArrayUp (m, size);
  return;
}


/*
 *                            A     B     C  D
 *                            \ \  / \  / /  |
 *                             \ /\  / \ /   |
 *                             /\ / \  / \   |
 *                            / /\   /\    \ |
 *                           //   \ /   \   \|
 *                          E      F     G   H
 *
 *    E  F  G  H
 * A  0  1  1  0    In this example, parent 'A' has children 'F' and 'G',
 * B  1  0  0  1                     parent 'B' has children 'E' and 'H',
 * C  1  1  0  0                     parent 'C' has children 'E' and 'F',
 * D  0  0  0  1                 and parent 'D' has child 'H'.
 *
 * "row barycenter" is average of the positions of a node's children.
 * "column barycenter" is average of the positions of a node's parents.
 *
 * pos:     0    1    2    3 ->
 *                              Row
 *       |  E    F    G    H  | Bary-y
 * ------+--------------------+-----
 * A     |  0    1    1    0  | 1.5  ==> 1+2/2 = 1.5
 * B     |  1    0    0    1  | 1.5  ==>   3/2 = 1.5
 * C     |  1    1    0    0  | 0.5  ==> 0+1/2 = 0.5
 * D     |  0    0    0    1  | 3.0  ==>   3/1 = 3
 * ------+--------------------+-----
 * Col   | 1.5  1.0  0.0  2.0 |
 * Baryx |                    |
 */


/* crossing edge count algorithm */

static inline int matrixCrossings (matrix * m);
static inline int matrixCrossingsInRows (matrix * m, int j, int k);

/* matrix and number of levels */
static inline unsigned int
matrixArrayCrossings (matrix * array[], int size)
{
  int i = 0;
  unsigned int sum = 0;
  unsigned int sumi = 0;
  /* for every level */
  sum = 0;
  for (i = 0; i < size; i++)
    {
      sumi = matrixCrossings (array[i]);
      sum = sum + sumi;
      if (lmgdebug)
	{
	  printf ("%s(): level %d has %d crossings total is %d\n",
		  __FUNCTION__, i, sumi, sum);
	}
      /* -1 is unknown number of crossings */
      if (current_edge_crossings != UINT_MAX)
	{
	  /* only need to know if layer has more crossings */
	  if (sum > current_edge_crossings)
	    {
	      break;
	    }
	}
    }
  return sum;

}

/* return the number of crossing edges between the nodes represented in
   the matrix m */
static inline int
matrixCrossings (matrix * m)
{
  int j = 0;
  int k = 0;
  int p = 0;
  int sum = 0;
  int sumi = 0;
  if (lmgdebug)
    {
      printf ("%s(): matrix has %d rows %d columns\n", __FUNCTION__, m->rows,
	      m->columns);
    }
  sum = 0;
  p = m->rows;
  for (j = 0; j < p - 1; j++)
    {
      for (k = j + 1; k < p; k++)
	{
	  sumi = matrixCrossingsInRows (m, j, k);
	  sum = sum + sumi;
	  if (lmgdebug && 0)
	    {
	      printf
		("%s(): crossings in rows is %d total %d for %d and %d\n",
		 __FUNCTION__, sumi, sum, j, k);
	    }
	}
    }
  return sum;
}

/* return the number of crossing edges between the nodes represented in
   rows j and k */
static inline int
matrixCrossingsInRows (matrix * m, int j, int k)
{
  int a;
  int b;
  int q;
  int sumi = 0;
  int sum = 0;
  sum = 0;
  q = m->columns;
  for (a = 0; a < q - 1; a++)
    {
      for (b = a + 1; b < q; b++)
	{
	  /* element at row,column 1/0 */
	  sumi = matrixAt (m, j, b) * matrixAt (m, k, a);
	  sum = sum + sumi;
	}
    }
  return sum;
}

/* misc */

/* crate a matrix with the given dimension (# of rows and columns).
   elements are initialized with integer zero */
static matrix *
matrixCreate (int r, int c)
{
  matrix *res;
  int i;
  int nbyte = 9;
  res = lmg_malloc (sizeof (matrix));
  res->rows = r;
  res->columns = c;
  res->logical_to_physical_row_number =
    lmg_malloc (r * sizeof (redirection_element));
  res->logical_to_physical_col_number =
    lmg_malloc (c * sizeof (redirection_element));

  /* edge bits */
  nbyte = ((r * c) + CHAR_BIT) / CHAR_BIT;
  res->elembits = lmg_malloc (nbyte * sizeof (unsigned char));

  for (i = 0; i < r; i++)
    {
      res->logical_to_physical_row_number[i].redirection = i;
      res->logical_to_physical_row_number[i].n = NULL;
      res->logical_to_physical_row_number[i].barycenter = 0.0;
    }
  for (i = 0; i < c; i++)
    {
      res->logical_to_physical_col_number[i].redirection = i;
      res->logical_to_physical_col_number[i].n = NULL;
      res->logical_to_physical_col_number[i].barycenter = 0.0;
    }
  return res;
}

/* End */
